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ABSTRACT 

In  this  paper,  we  propose  a  new  approach  to  wavelet-based  de- 
convolution.  Roughly  speaking,  the  algorithm  comprises  Fourier- 
domain  system  inversion  followed  by  wavelet-domain  noise  sup¬ 
pression.  Our  approach  subsumes  a  number  of  other  wavelet-based 
deconvolution  methods.  In  contrast  to  other  wavelet-based  ap¬ 
proaches,  however,  we  employ  a  regularized  inverse  filter,  which 
allows  the  algorithm  to  operate  even  when  the  inverse  system  is  ill- 
conditioned  or  non-invertible.  Using  a  mean-square-error  metric, 
we  strike  an  optimal  balance  between  Fourier-domain  and  wavelet- 
domain  regularization.  The  result  is  a  fast  deconvolution  algorithm 
ideally  suited  to  signals  and  images  with  edges  and  other  singu¬ 
larities.  In  simulations  with  real  data,  the  algorithm  outperforms 
the  LTI  Wiener  filter  and  other  wavelet-based  deconvolution  algo¬ 
rithms  in  terms  of  both  visual  quality  and  MSE  performance. 

1.  INTRODUCTION 

Deconvolution  is  a  recurring  theme  in  a  wide  variety  of  signal  and 
image  processing  problems,  from  channel  equalization  [1]  to  im¬ 
age  restoration  [2],  Often,  the  distortion  introduced  by  a  measure¬ 
ment  device  can  be  modeled  as  a  convolution  of  the  desired  data 
with  the  impulse  response  of  the  device.  Deconvolution,  then,  cor¬ 
responds  to  inverting  the  effects  of  the  distortions.  Unfortunately, 
the  measured  signal  is  usually  also  corrupted  by  noise,  which  com¬ 
plicates  the  process  of  deconvolution. 

In  its  simplest  form,  the  1-d  deconvolution  problem  runs  as 
follows.  The  desired  signal  x  is  input  to  a  known  linear  time- 
invariant  (LTI)  system  having  impulse  response  h.  White  Gaus¬ 
sian  noise  n  of  variance  a2  corrupts  the  output  of  the  system.  We 
measure  the  result  y(t)  :=  ( h  *  x)(t)  +  n(t).  In  the  Fourier  do¬ 
main,  we  have  Y(f)  =  H(f)  X(f)  +  N(f).  Given  y,  we  seek  to 
estimate  x. 

If  the  system  frequency  response  H(f)  has  no  zeros, 
then  we  can  obtain  an  unbiased  estimate  of  x  as  X(f)  := 
=  X{f)  +  H-\f)N(f).  However,  if  H(f)  be¬ 
comes  small  at  any  frequency,  then  enormous  noise  amplification 
resuits,  yielding  an  infinite-variance,  useless  estimate. 

In  situations  involving  such  ill-posed  systems,  some  amount  of 
regularization  becomes  essential.  Regularization  reduces  the  vari¬ 
ance  of  the  signal  estimate  (noise  reduction)  in  exchange  for  an  in¬ 
crease  in  bias  (signal  distortion).  When  x  is  wide-sense  stationary 
(WSS),  the  LTI  Wiener  filter  provides  the  optimal  regularization  in 
the  minimum  mean-squared-error  (MSE)  sense  [2]. 

Unfortunately,  the  signals  and  images  appearing  in  many  im¬ 
portant  applications  contain  information-bearing  edges  and  ridges. 
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The  LTI  Wiener  filter  is  inappropriate  for  such  non-stationary  sig¬ 
nals,  for  since  it  reduces  noise  by  smoothing  all  signal  components 
uniformly;  it  can  smear  local  features  such  as  edges.  The  root  of 
this  problem  lies  in  the  fact  that  while  the  underlying  Fourier  basis 
of  the  LTI  Wiener  hlter  matches  an  LTI  system  h ,  it  does  not  match 
a  non-stationary  signal  x. 

Wavelets,  on  the  other  hand,  provide  a  basis  matched  to  a 
large  class  of  non-stationary  signals  [3,4],  Unlike  the  strict  fre¬ 
quency  localization  of  the  Fourier  basis,  a  wavelet  basis  localizes 
in  both  time  and  frequency,  and  so  can  effectively  track  signal 
non-stationarities.  This  matching  property  has  been  leveraged  into 
powerful  algorithms  for  noise  reduction  that  simply  threshold  the 
wavelet  representation  of  the  noisy  signal  [5, 6].  Wavelet  denois- 
ing  is  a  spatially  adaptive  processing  (it  smoothes  more  in  smooth 
regions  of  the  signal)  ideally  suited  to  signals  with  edges  and  other 
singularities. 

The  fact  that  LTI  systems  are  matched  by  one  basis  (Fourier) 
but  non-stationary  signals  by  another  (wavelet)  inspires  a  hybrid 
approach  to  deconvolution:  (i)  invert  the  convolution  operation  in 
the  Fourier  domain  and  then  (ii)  regularize  (denoise)  in  the  wavelet 
domain.  Such  an  approach  has  been  followed  by  Donoho  [7], 
Nowak  [8],  and  Mallat  [3,  pp.  456-461]  with  considerable  success. 

However,  current  wavelet-based  deconvolution  schemes  can¬ 
not  deal  with  ill-conditioned  systems  h,  since  they  entrust  all  of 
the  regularization  to  their  wavelet  denoising  post-processing  step. 
Systems  with  zeroes  in  the  frequency  response  will  present  noise 
of  infinite  variance  to  the  denoising  step,  destroying  the  signal  es¬ 
timate. 

In  this  paper,  we  propose  an  improved  hybrid  Fourier/wavelet 
deconvolution  algorithm  suitable  for  use  with  ill-conditioned  sys¬ 
tems.  The  basic  idea  is  simple:  employ  both  Fourier-domain 
(Wiener-like)  and  wavelet-domain  regularization.  With  this  tan¬ 
dem  processing,  we  can  keep  the  Fourier-domain  regularization 
(and  its  corresponding  smearing  distortions)  to  the  minimum  re¬ 
quired  to  make  the  system  transfer  function  well-posed;  the  bulk 
of  the  noise  removal  comes  in  the  wavelet  denoising  stage.  Using 
the  MSE  metric,  we  will  strike  an  optimal  balance  between  global 
and  local  processing.  Interestingly,  one  extreme  of  the  balance  is 
to  perform  no  Fourier-domain  regularization,  and  this  coincides  to 
the  approaches  of  [3,7,8].  Extension  to  multi-dimensional  data  is 
trivial. 

After  discussing  regularization  in  more  depth  in  Section  2  and 
previous  Fourier/wavelet  deconvolution  approaches  in  Section  3, 
we  present  our  improved  scheme  in  Section  4.  Illustrative  exam¬ 
ples  lie  in  Section  5.  We  close  with  conclusions  in  Section  6. 

2.  REGULARIZED  INVERSE  FILTERS 

Consider  a  zero-mean,  WSS  signal  x  with  power  spectral  density 
(PSD)  Px(f).  Given  the  general  deconvolution  problem  from  the 
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Introduction,  a  general  form  for  a  Fourier-domain-regularized  sig¬ 
nal  estimate  is  given  by  [9] 

Xa(f)  ■■=  Ga(f)Y(f)  (1) 


with 

GAf) 


1  ^  (  \H(f)\2Px(f)  \ 

H(f)J  \\H{f)\*Px(f )  + 


(2) 


The  regularization  parameter  a  controls  the  tradeoff  between  the 
amount  of  noise  suppression  and  the  amount  of  signal  distortion. 
Setting  a  =  0  gives  an  unbiased  but  noisy  estimate.  Setting 
a  =  o o  completely  suppresses  the  noise,  but  also  totally  distorts 
the  signal  (x<x>  =  0).  For  a  =  1,  (2)  corresponds  to  the  LTI 
Wiener  filter,  which  is  optimal  in  the  MSE  sense  for  a  Gaussian 
input  signal  x. 

Since  the  Fourier  basis  functions  underlying  any  LTI  filter 
have  spatial  support  over  the  entire  signal,  Ga  will  tend  to  smear 
non-stationarities  in  the  desired  signal,  such  as  edges  and  ridges. 
While  for  non-stationary  signals  we  could  just  solve  the  more  gen¬ 
eral  MSE  signal  estimation  problem  (time-varying  Wiener  filter¬ 
ing),  such  an  approach  would  have  no  special  structure  and  further¬ 
more  would  require  precise  knowledge  of  the  non-stationarities  in 
x. 


3.  WAVELETS  AND  DECONVOLUTION 

The  joint  time-frequency  analysis  of  the  wavelet  basis  efficiently 
captures  non-stationary  signal  features.  The  discrete  wavelet  trans¬ 
form  (DWT)  represents  a  1-d  signal  z  in  terms  of  shifted  versions 
of  a  low-pass  scaling  function  cp  and  shifted  and  dilated  versions 
of  a  prototype  bandpass  wavelet  function  ip  [3,4],  For  special 
choices  of  <p  and  ip,  the  functions  ipj,h(t)  2 ~3^2  ip(2~3t  —  It), 
<pj,k(t )  :=  2~3  <p(2~3t  —  fc),  j,k  £  Z  form  an  orthonomial  ba¬ 
sis,  and  we  have  the  representation  [3,4] 

_  jo  _ 

At)  =  E  jo, k  tfijojk  a)  +  E  E  Wj,hipj,k(t),  (3) 

k  j=  —  ook 

with  Uj,k  f  z(t)  <pj,k(t)  dt  and  Wj,k  f  z(t)  ipj,k(t)  dt.  For 

brevity,  we  will  collectively  refer  to  the  set  of  scaling  and  wavelet 
coefficients  as  { @ y , fc }  :=  {uj0,k,  wy*.}.  Multidimensional  DWTs 
are  easily  obtained  by  alternately  wavelet-transforming  along  each 
dimension  [3,4], 

The  DWT  enjoys  an  enviable  energy  compaction  property:  the 
energy  of  many  real-world  signals  compacts  into  just  a  few  large 
wavelet  coefficients,  while  white  noise  remains  disbursed  over  a 
large  number  of  small  coefficients.  This  disparity  can  be  exploited 
to  distinguish  signal  from  noise  and  has  given  rise  to  a  number  of 
powerful  denoising  techniques  based  on  simple  thresholding  [4— 
6]  that  can  suppress  noise  while  preserving  time-localized  signal 
structures. 

Wavelet  denoising  figures  prominently  in  a  number  of  recent 
advanced  deconvolution  algorithms  [3,  7,  8].  All  three  methods 
have  the  same  two  basic  steps  in  common: 

Inversion:  Compute  the  noisy  estimate  Xo  = 

This  inversion  necessarily  amplifies  noise  components  at 
frequencies  where  H(f)  is  small. 

Regularization  by  Wavelet  Denoising:  Compute  the  DWT  of 
xo,  and  then  threshold  and  invert  the  DWT  to  obtain  the 
final  signal  estimate  x.  Note  that  the  Inversion  step  colors 
the  white  corrupting  noise  n;  hence  scale-dependent  thresh¬ 
olds  [6]  should  be  employed. 


Regularized 

_ 

Wavelet 

Inverse  Ga 

^  $Oi  ^ 

Denoising 

Figure  1:  Wavelet-based  regularized  deconvolution  (WaRD):  par¬ 
tially  regularized  inverse  bitering  following  by  wavelet  denoising. 


4.  IMPROVED  WAVELET-BASED  DECONVOLUTION 
WITH  REGULARIZED  INVERSE 

The  current  Fourier/wavelet  deconvolution  approaches  of  [3, 7, 8] 
completely  decouple  the  inversion  and  regularization  processes. 
Unfortunately,  when  the  system  h  is  very  ill-conditioned  (or  sim¬ 
ply  non-invertible),  any  attempt  at  inversion  will  amplify  the  cor¬ 
rupting  noise  to  an  extent  that  it  will  obliterate  the  desired  signal. 
No  amount  of  wavelet  denoising  can  rescue  us  in  this  case. 

Note,  however,  the  sensitivity  of  the  inversion  process  to  regu¬ 
larization.  A  minute  amount  of  regularization  (small  a  in  (2))  can 
lead  to  a  huge  reduction  in  the  degree  of  noise  amplification  -  all 
this  at  the  expense  of  only  a  slight  increase  in  the  signal  distortion. 
This  realization  motivates  us  to  replace  the  Inversion  step  of  the  al¬ 
gorithms  of  [3,7, 8]  with  a  Regularized  Inversion  step  (see  Figure 
l).1  We  call  the  resulting  algorithm  wavelet-domain  regularized 
deconvolution  (WaRD). 

But  how  to  pick  the  right  value  for  the  regularization  param¬ 
eter  a?  The  tradeoff  is  clear:  On  one  hand,  since  regularization 
smears  non-stationary  signal  features  like  edges  and  ridges,  we 
would  prefer  a  as  small  as  possible.  On  the  other  hand,  large  a 
prevents  excessive  noise  amplification  during  inversion  which  aids 
the  wavelet  denoising. 

To  be  more  precise,  we  will  determine  the  optimal  regular¬ 
ization  parameter  for  the  WaRD  system  by  minimizing  the  over¬ 
all  MSE.  The  deconvolution  MSE  consists  of  the  signal  distortion 
due  to  Fourier-domain  regularization  and  the  error  due  to  wavelet- 
domain  denoising: 

MSE(a)  =  f[l-Ga  (/)  H(f)]2  Px  (/)  df 

+  E  min(|flj.*(&»)|ai  °j(Q))  ■  (4) 

j,k 

Here  6j,k(xa )  denotes  the  wavelet  coefficients  of  xa,  ct|(q)  is 
the  variance  of  the  wavelet-domain  noise  at  scale  j,  and  Px  (/)  = 

l*(/)|2- 

The  first  term  in  MSE(ct)  is  an  estimate  of  the  distortion  in 
the  input  signal  due  to  the  regularized  Fourier-domain  inverse  [9]. 
This  distortion  is  an  increasing  function  of  a.  The  second  term  is 
an  estimate  of  the  error  due  to  ideal  wavelet  domain  hard  thresh¬ 
olding  [5],  Ideal  thresholding  consists  of  keeping  a  noisy  wavelet 
coefficient  only  if  the  signal  power  in  that  coefficient  is  greater 
than  the  noise  power.  Otherwise,  the  coefficient  is  set  to  zero. 
(Ideal  thresholding  assumes  that  the  signals  under  consideration 
are  known.)  This  error  is  a  decreasing  function  of  a.  The  optimal 
regularization  parameter,  denoted  by  a*,  corresponds  to  the  min¬ 
imum  of  MSE  (a).  Note  that  a*  depends  on  the  signal,  system, 
and  noise  power. 

The  existing  Fourier/wavelet  deconvolution  algorithms  of  [3, 
7,  8]  can  be  interpreted  as  special  cases  of  WaRD  with  a  —  0. 

1Note  that  the  Fourier-domain-regularized  inverse  (2)  was  derived  for 
WSS  signals  x  only.  With  non-stationary  x.  we  replace  Px  by  the  time- 
averaged  spectrum  of  x.  In  practice,  we  set  Px  (f)  =  X (/)  | 2 . 


However,  as  mentioned  earlier,  these  methods  are  in  general  not 
applicable  when  h  is  not  invertible.  Even  when  h  is  invertible, 
WaRD  will  outperform  these  methods,  since  the  value  a  =  0  is 
included  in  the  search-space  for  the  optimal  a*. 

The  computational  cost  of  deconvolving  an  M-point  signal 
will  be  dominated  by  the  0(M  log  M)  cost  of  the  FFT  inverse- 
filter  implementation.  The  second,  wavelet  denoising  step  con¬ 
sumes  only  O(M)  computations. 

The  wavelet  denoising  step  of  the  WaRD  algorithm  can  be  ex¬ 
tended  in  several  ways.  First,  since  the  standard  DWT  is  not  shift- 
invariant,  shifts  of  y  will  result  in  different  estimates  x.  Employing 
a  redundant,  shift-invariant  DWT  will  both  yield  a  shift-invariant 
algorithm  as  well  as  improve  the  denoising  performance  substan¬ 
tially  [4],  all  at  no  significant  increase  in  the  overall  computational 
cost.  The  recently  proposed  complex  and  “almost  shift-invariant” 
DWT  of  [10]  would  yield  similar  results  at  a  reduced  computa¬ 
tional  cost.  Finally,  instead  of  a  threshold,  we  can  apply  a  Wiener 
filter  to  the  wavelet  coefficients  [11,  12], 2  Such  processing  has 
been  shown  to  outperform  simple  thresholding  for  denoising  finite 
samples  of  data. 

Finally,  note  that  WaRD  extends  trivially  to  higher  dimensions 
using  the  appropriate  Fourier  and  wavelet  transforms. 


dashed  lines  in  Figure  2(b).  The  denoising  step  was  implemented 
by  hard-thresholding  the  coefficients  of  a  shift-invariant  DWT  us¬ 
ing  the  best  wavelet  packet  basis.  This  algorithm  outperforms  the 
methods  of  [7, 8]  typically,  as  well  as  the  standard  Wiener  filter  in 
this  case  (see  Figure  2(e)). 

Figure  2(f)  plots  the  WaRD  obtained  using  a*  =  0.06. 
Wavelet-domain  Wiener  filtering  was  applied  to  the  coefficients 
of  a  shift-invariant  DWT  for  the  denoising  stage.  The  WaRD  out¬ 
performs  the  other  algorithms  in  terms  of  both  visual  quality  and 
MSE  performance. 

Next,  we  consider  image  restoration  using  WaRD  (Re WaRD). 
The  input  x  is  the  256  x  256  Lena  image  (normalized  to  zero  mean 
and  unit  energy)  and  the  discrete-time  system  response  h  is  a  2-d, 
4-point  smoother  [111  1]T  [1  1  1  1].  Such  a  response  is 
commonly  used  as  a  model  for  blurring  due  to  a  square  scanning 
aperture  such  as  in  a  CCD  camera  [2],  The  noise  variance  was 
set  to  a 2 * *  =  4  x  1CP7.  Figure  3  illustrates  the  desired  x,  the  ob¬ 
served  y,  the  Wiener  filter  estimate  xi,  and  the  WaRD  estimate  for 
a*  =  0.27.  The  methods  of  [3,7,8]  are  not  applicable  in  this  situ¬ 
ation,  due  to  the  many  zeros  in  H(fx,  fy).  (The  Wiener  filter  out¬ 
performed  the  wavelet  packet  method  in  this  case.)  WaRD  clearly 
outperforms  Wiener  filtering  in  both  visual  quality  and  MSE. 


5.  EXAMPLES 


6.  CONCLUSIONS 


To  illustrate  the  performance  of  the  WaRD  algorithm,  we  will  per¬ 
form  simulations  in  1-d  and  2-d. 

We  first  compare  the  WaRD  with  other  methods  for  the  1-d 
deconvolution  problem  presented  in  the  Introduction.  In  order  to 
observe  the  behavior  of  each  method  for  both  smooth  and  edgy 
regions,  we  take  for  x  a  concatenation  of  Donoho’s  Blocks  and 
Heavisine  signals  [13]  (normalized  to  be  zero  mean  and  unit  en¬ 
ergy).  We  employ  Daubechies  length-8  wavelets  throughout.  Fig¬ 
ure  2(a)  depicts  the  signal.  For  the  system,  we  take  the  example 
of  [3,  pp.  459] 


1, 

2-4|/|, 


I/I  €  [0,0.25] 
|/|  e  (0.25,0.5] 


(5) 


with  frequency  /  normalized  to  (—0.5,  0.5]  (see  Figure  2(b)).  The 
corrupting  noise  variance  was  a1  =  4  x  10~6.  Figure  2(c)  plots 
the  blurred,  noisy  signal. 

The  Wiener  filter  estimate  (Figure  2(d))  was  implemented  us¬ 
ing  ct  =  1  andPa,(/)  =  |X(/)|2  in  (2).  The  Wiener  filter  bases  its 
deconvolution  on  the  signal-to-noise  ratio  at  each  frequency.  How¬ 
ever.  this  is  inappropriate  for  non-stationary  signals,  since  their 
frequency  content  changes  with  time.  The  deconvolution  suffers 
in  response. 

The  methods  of  [7,  8]  fail  in  this  case,  due  to  the  null  in  the 
frequency  response  of  the  system  at  77(0.5)  =  0. 

Since  H(f)  decays  to  zero  at  /  =  0.5  so  slowly,  the 
wavelet  packet  decotwolution  method  of  [3,  pp.  458-461]  remains 
applicable.3  This  method  adapts  a  wavelet  packet  basis  to  the  col¬ 
ored  noise  in  the  inverted  data  X (/)  =  Y (/).  The  fre¬ 

quency  splits  of  the  wavelet  packet  basis  for  H  are  shown  with 

2  Do  not  confuse  wavelet-domain  Wiener  filtering  with  the  Fourier- 
domain  Wiener  filtering  discussed  above. 

3Even  though  H(f )  is  not  invertible,  the  wavelet  packet  deconvolution 

method  does  not  fail,  because  the  location  and  the  order  of  the  zero  of 

H  are  such  that  only  a  few,  high-frequency  signal  wavelet  coefficients  are 

obliterated  by  the  infinite  noise  amplification  at  /  =  0.5.  This  method  will 

fail  when  11(f)  has  zeros  of  arbitrary  location  and  order,  however. 


In  this  paper,  we  have  proposed  an  efficient  multi-scale  deconvo¬ 
lution  algorithm  that  optimally  combines  Fourier-domain  regular¬ 
ized  inversion  and  wavelet-domain  denoising.  For  non-stationary 
signals,  the  WaRD  outperforms  the  LTI  Wiener  filter  and  other 
wavelet-based  deconvolution  algorithms  in  terms  of  both  visual 
quality  and  MSE  performance.  Furthermore,  it  continues  to  pro¬ 
vide  a  good  estimate  of  the  original  signal  even  when  the  system 
response  is  ill-conditioned.  All  of  this  in  an  algorithm  of  compu¬ 
tational  complexity  no  greater  than  an  FFT. 

For  a  given  problem  setup,  the  optimal  value  of  the  regular¬ 
ization  parameter  a  depends  on  the  signal,  the  system  response, 
and  the  noise  level.  Fortunately,  final  performance  was  observed 
to  be  quite  insensitive  to  the  exact  value,  as  long  as  we  choose  al¬ 
pha  sufficiently  positive.  As  a  guide,  in  simulations  spanning  many 
real-world  images  and  convolution  systems,  a*  was  almost  always 
lay  in  the  range  [0.2,  0.3],  In  the  1-d  example  above,  the  optimal 
a  turned  out  to  be  small  (a*  =  0.06)  because  the  test  signal  com¬ 
pacted  very  well  in  the  wavelet  domain.  However,  choosing  a  a* 
in  the  range  [0.2,  0.3]  gave  near-optimal  results. 

There  are  several  avenues  for  future  WaRD  related  research. 
We  are  developing  methods  to  estimate  the  WaRD  MSE,  and  thus 
the  optimal  a*,  without  prior  knowledge  of  the  signal  and  noise 
power.  Further,  we  are  investigating  the  gains  possible  using  a 
"best  basis”  adapted  to  the  signal  instead  of  the  noise  as  suggested 
in  [3,  pp.  458-^-61]. 
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